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Abstract: A prototype silicon-tungsten electromagnetic calorimeter for an ILC detector was 
tested in 2007 at the CERN SPS test beam. Data were collected with electron and hadron beams 
in the energy range 8 to 80 GeV. The analysis described here focuses on the interactions of pions 
in the calorimeter. One of the main objectives of the CALICE program is to validate the Monte 
Carlo tools available for the design of a full-sized detector. The interactions of pions in the Si-W 
calorimeter are therefore confronted with the predictions of various physical models implemented 
in the GEANT4 simulation framework. 

KEYWORDS: Detector physics: concepts, processes, methods, modelling and simulations; 
Calorimeter methods; Detector modelling and simulations I (interaction of radiation with matter, 
interaction of photons with matter, interaction of hadrons with matter, etc.) . 
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1. Introduction 

The current designs of the calorimetry for detectors at the ILC are in large part driven by the 
demands of jet energy reconstruction. A key aim is to achieve a relative energy resolution of 
~ 30% / -y/Zs/GeV ; this precision would for example permit the reconstruction of the hadronic 
decays of the W and Z bosons with a precision comparable with their natural widths, and would 
thus allow the W and Z to be distinguished in their hadronic (two-jet) decay modes. 

This target for jet energy resolution is roughly a factor two better than achieved in previ- 
ous detectors. The most promising way to achieve it is believed to be through the "particle 
flow" approach Jj]]. The idea is that particles of different types in jets should be reconstructed in 
the different parts of the detector where they can be measured most precisely: the charged particles 
in the tracking system, photons in the electromagnetic calorimeter (ECAL) and neutral hadrons in 
a combination of the ECAL and the hadronic calorimeter (HCAL). The key to this concept is to 
minimise confusion in the pattern recognition in the calorimeters, and to achieve this, high spatial 
granularity is required. The optimum design of the ILC detectors can be addressed by Monte Carlo 
simulation, but in order to do this, it is crucial first to validate the Monte Carlo tools against data. 

A first round of beam tests was performed at DESY and CERN in summer 2006, followed by 
more complete tests in 2007. These tests were performed with a combined system of a silicon- 
tungsten (Si-W) ECAL, followed by a hadronic calorimeter and then a coarser tail catcher, the 
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latter two both built from a scintillator-iron sandwich structure. Since the ECAL has a length 
~ 1 .OAim. , it will seldom fully contain a hadronic shower. However, more than half of hadronic 
showers will start their development in the ECAL, so it is important to investigate the properties 
of hadronic showers initiated in the (mainly tungsten) material of the ECAL. The ECAL also has 
some advantages for the study of hadronic shower properties. For example, the interaction point 
can be identified with good precision. The small ratio of Xq / Ant. (a factor ~ 3 smaller for tungsten 
than for iron) means that the electromagnetic products of the first interaction can be expected to 
shower rapidly after the interaction point, and we can hope to exploit the fine granularity of the 
detector to separate them from the hadronic products. This means that the longitudinal shower 
development after the interaction point is potentially a particularly interesting quantity. 

In Sect. H we briefly describe the main features of the Si-W calorimeter and in Sect. [| we 
outline the test beam setup at CERN. Then in Sect. ^| we review the Monte Carlo models with 
which the data are confronted. In Sect. || we explain the selection of data for comparison with 
simulation in Sect. ||. 

2. The CALICE Si-W ECAL prototype 

The ECAL prototype used in this study is a silicon-tungsten sampling calorimeter, made of 30 
readout layers [^] . The active detectors were silicon wafers segmented into a 6 x 6 array of diode 
pads, each with a size of 1 x 1 cm 2 . The mechanical structure was built from tungsten sheets 
wrapped in carbon fibre. Between each sheet of tungsten, the structure contained alveolar slots, 
into which detector slabs were inserted, each consisting of a further tungsten sheet sandwiched 
between two layers of sensors. The prototype was constructed from three stacks, each composed of 
ten layers of alternating tungsten and silicon, and each stack having a different tungsten thickness: 
1.4 mm or 0.4X per layer in the first stack, 2.8 mm or 0.8X per layer in the second stack and 
4.2 mm or \.2Xq per layer in the rear stack. In terms of interaction lengths, these tungsten layers 
each contribute ~ 1.4% Aim., ~ 2.7% Aim. an d ~ 4.1% Ant. respectively. Other materials (mainly 
the carbon-fibre-epoxy mechanical structure, PCBs and the silicon itself) contribute alternately 
~ 0.2%Ai nt . and ~ 1.0%Ai nt . between successive samplings. The overall thickness is about 20 cm, 
corresponding to ~ 24.6Xq or ~ l.OAm. at normal incidence. 

A detailed description of the prototype's hardware and of its commissioning in test beams at 
DESY and CERN in 2006, can be found in ref. to which the reader is referred for much fuller 
details. 

The full prototype consists of a 3x3 array of wafers (i.e. an 18x18 array of 1 cm 2 pads) in 
each layer. The version of the prototype tested at CERN in 2007 consisted of 30 layers, of which 
initially the first twelve were instrumented with a 3x2 array of wafers, and the remaining eighteen 
were fully equipped with a 3x3 array. Later in the 2007 run six more layers were completed, 
leaving only the first six in a 3x2 configuration. 

In offline analysis, the raw hit energies in each cell have their pedestals subtracted, and are 
converted from raw ADC counts into Minimum Ionising Particle equivalents (MIPs), as explained 
in ref. [g]. This gain correction is derived by finding the most probable energy deposition in each 
cell in data recorded with a high energy muon beam with a typical relative precision of ~ 0.5% 
per pad. 



-2- 



TOP 



Cerenkov 
11000 



Scl DC 3 Sc2 Veto Sc3 DC2 
15 29 



123 I 407 



DC1 

58 



ECAL HCAL TCatcher 
201 1223 1458 



-2417.5 
-2572 - g47 



-29 



1587.5 2327.5 4018 4812 



FRONT Veto is 1000x1000, with a hole of 200x200 mm 
Scl andSc3 are 100x100 
Sc2 is 200x200 



Mcl is 1000x1000 All distances are in mm 



Figure 1. Schematic layout for the CALICE beam tests at CERN in 2007 (not to scale). Dimensions are 
indicated in mm. 



3. The test beam 

In this paper we report on some of the data taken in 2007 in the CERN H6 test beam [g]. The 
layout of the CALICE calorimeters in the test beam is shown schematically in figure [|. The three 
calorimeters used were the Si-W ECAL, followed by an HCAL prototype [Q] (using iron as the 
absorber and 38 layers of scintillator tiles with analogue readout as the active medium, using tile 
size 3x3 cm in the shower core) and a tail catcher (TCMT - also an iron calorimeter with sixteen 
layers of 5 cm wide scintillator strips). The thickness of the iron sheets was ~ 18 mm in the HCAL, 
~ 20 mm for the first eight layers of the TCMT, and ~ 100 mm in the rear section of the TCMT. 
The ECAL and HCAL were mounted on a movable stage, providing the possibility to translate and 
rotate the calorimeters with respect to the beam. 

Upstream of the calorimeters were scintillation counters used for triggering. A muon veto 
counter was also installed downstream of the TCMT. The standard trigger used for the present 
analysis required a coincidence between any two of the upstream counters Scl, Sc2 and Sc3. In 
addition, three sets of wire chambers were operated to permit the beam position and direction to be 
measured upstream of the calorimeters. A Cerenkov counter was used in threshold mode to assist 
particle identification. 

Data were recorded in 2007 using electron, hadron and muon beams with energies in the range 
8 to 180 GeV. A variety of different calorimeter translations and rotations were used, and in all 
~200 million triggers were recorded, including calibration data. For the present study, we used 
data corresponding to a configuration with the beam impinging at normal incidence close to the 
centre of the calorimeters. The analysis uses %~ beams at energies of 8, 10, 12, 15 and 20 GeV, 
and 7T + runs at 30, 50 and 80 GeV. Each run comprised typically ~200,000 triggers and ~ 100,000 
hadronic events. For the %~ runs the pressure in the Cerenkov counter was set so that electrons 
could be vetoed, while for the n + runs, it was set so as to separate protons from pions. 

4. Monte Carlo simulation 

The main purpose of the present analysis is to confront the data with Monte Carlo simulations, 
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and specifically to test the different physical models available. These simulations were carried out 
in the framework of GEANT4 [S], using version 4.9.3, released in December 2009. The results 
presented are based on samples of typically 50000 simulated events. 

In order to describe the geometry of the CALICE prototypes within GEANT4 we used the pro- 
gram Mokka [§], which is also capable of simulating full ILC detector geometries. For the current 
study, the simulation of the ECAL is the most important. Individual silicon sensors, with their pad 
structure and guard rings, are represented in Mokka, as well as the tungsten radiator, and other 
passive materials such as the carbon-fibre and epoxy of the support structure and the PCBs used for 
readout. The layer-by-layer staggering of the wafer positions is simulated. Similar representations 
of the HCAL and TCMT are also included. The detectors upstream of the calorimeter systems 
(Cerenkov, scintillators and tracking chambers) are also modelled by Mokka. 

During particle transport in GEANT4, Mokka records the ionisation energy deposited in sen- 
sitive detectors (Silicon pads, scintillators etc.). A simple simulation of the noise contribution in 
each cell is performed, as outlined in ref. [{7|]. At the end of the simulation, the summed energy 
deposits in each cell are finally converted into MIPs. The conversion factor is based on matching 
the minimum ionising energy peak in data and simulation. This procedure is estimated to lead to 
a ~ 1% systematic uncertainty in the relative energy scales of the data and simulation. As a first 
step in the analysis, a hit energy cut of 0.6 MIPs is imposed in both simulation and data, in order 
to eliminate most noise-only hits. The hits are also required to occur within 150 ns of the beam 
arrival time at the trigger counters, in order to emulate the behaviour of the readout electronics. 

In the simulation, pions are simulated starting with a Gaussian transverse profile in x and y, x 
at a point ~60 m upstream of the calorimeters, so that they pass through, and potentially interact 
in, the Cerenkov counter, trigger scintillators, tracking chambers and intervening air. Hits are 
recorded in the scintillators and tracking chambers, and energies deposited in the scintillators are 
used as the basis of a simple simulation of the trigger. The beam profile in the ECAL, defined by the 
distributions of the shower barycentre in x and y, is then compared between data and simulation. 
The input Gaussian parameters are adjusted until satisfactory agreement is achieved. A typical 
example of the comparison after tuning is shown in figure ||. 

GEANT4 provides the user with a number of "physics lists" — combinations of models se- 
lected to simulate particle interactions for different species at different energies. In our study we 
are particularly concerned with the hadronic interactions for pions. Typically these physics lists use 
different models in different energy ranges, with smooth interpolation between models achieved by 
allowing the energy ranges to overlap, with smoothly varying random selection between the models 
in the changeover region. Detailed descriptions of the models can be found in ref. [JsJ] . The fol- 
lowing physics lists have been chosen for our investigation, based on advice and recommendations 
from the GEANT4 authors: 

LHEP Uses the LEP (low energy parametrised) and HEP (high energy parametrised) models, mak- 
ing a transition between the two over the range 25-55 GeV. These are essentially a recoding 
into C++ of the GHEISHA model [Q] extensively used in earlier simulations using GEANT3, 

'The CALICE coordinate system is defined with z along the nominal beam direction, y vertical, and x horizontal. 
The origin is at the centre of the exit face of the most downstream tracking chamber. 
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Figure 2. Example of the beam profiles in x andy (based on the shower centroid reconstructed in the ECAL) 
observed with a 20 GeV 7t~ beam. Data (points with error bars) are compared with the tuned simulation 
(solid histogram). The distributions are normalised to the same numbers of events. 



for example. It is still widely used, even though it is not regarded as the state-of-the-art 
choice. 



FTFP_BERT Uses the GEANT4 implementation of the Bertini cascade model 2 [|0|] for low ener- 



gies, making a transition to the FTFP model, based on the GEANT4 implementation of the 
Fritiof diffractive string model [pT|], at pion energies between 4 and 5 GeV. We have also 
studied the FTFP_BERT_TRV list, which has the same model content as FTFP_BERT, but 
with a higher transition energy, in the region 6-8 GeV. In practice, for the data presented 
here, we find no significant differences between these two physics lists, so we only show 
FTFP_BERT here. 

QGSP_BERT Uses the Bertini model [|l^] at low energies, making a transition to the LEP (GHEI SHA) 
model between 9.5 and 9.9 GeV, and a further transition to the QGSP model between 12 and 
25 GeV. QGSP uses a GEANT4 implementation of a string model [ |12[ ] for the high energy 
interaction, supplemented by the GEANT4 precompound model[|l3|] describing deexcitation 
of the nucleus. We have also studied the QGSP_BERT_HP list, which has the same model 
content as QGSP_BERT, combined with a more accurate treatment of low energy neutron 
scattering. For the distributions presented in this paper we find no significant benefit in using 
the QGSP_BERT_HP list, and hence we only present QGSP_BERT. 

QGSP_BERT_TRV has the same model content as QGSP_BERT, but with the transition between 
LHEP and QGSP occurring at lower energy, in the region 10-15 GeV, in order to reduce the 
reliance on LHEP. 

QGSP_FTFP_BERT Similar to QGSP_BERT, with different transition energies, and using FTFP 
instead of LEP in the intermediate region. 



QGS_BIC Uses the GEANT4 binary cascade model (BIC) [14] at the lowest energies, then LEP 



in the intermediate region, and QGS at high energies. The BIC model is also used for the 



2 This model incorporates the Bertini intra-nuclear cascade model with excitons, a pre-equilibrium model, a nucleus 
explosion model, a fission model, and an evaporation model. 
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rescattering of secondaries in this case (denoted QGSB below) and the GEANT4 precom- 
pound model is used to describe deexcitation of the nucleus. 

QGSP_BIC In this case the BIC model is not used for pions (only for neutrons and protons), and 
there is no rescattering of secondaries. LEP is used in the low energy region, and QGSP at 
high energies. 

FTF_BIC Uses the BIC model at low energies, with a transition to Fritiof between 4 and 5 GeV. 
The BIC model is again used for the rescattering of secondaries in this case (denoted FTFB 
below) and the GEANT4 precompound model is used to describe deexcitation of the nucleus. 

The physics content of these models for pions is summarised in table [I[ The QGSP_BERT physics 
list is generally favoured by the LHC general purpose detectors for calorimetry applications [|l5|], 
based mainly on data from test beams of higher energies than those relevant for CALICE. In our 
analysis, we expect to be particularly sensitive to the model(s) used for the primary interaction at 
the energy of the incident beam. This is indicated for each physics list in table ||. 

Table 1. The table shows the physics models invoked for pion inelastic interactions in each physics list. 
Where ranges overlap, GEANT 4 chooses randomly between models, with probabilities varying linearly with 
energy over the range of overlap. 



Physics List 


Model content (for 


FTFP_BERT 

QGSP_BERT 

QGSP_BERT_TRV 

QGSP_FTFP_BERT 

QGS_BIC 

QGSP_BIC 

FTF_BIC 

LHEP 


Bertini (0-5 GeV); FTFP (>4 GeV) 

Bertini (0-9.9 GeV); LEP (9.5-25 GeV); QGSP (> 12 GeV) 
Bertini (0-9.9 GeV); LEP (9.5-15 GeV); QGSP (> 10 GeV) 
Bertini (0-8 GeV); FTFP (6-25 GeV); QGSP (>12GeV) 
BIC (0-1.3 GeV); LEP (1.2-25 GeV); QGSB (>12GeV) 
LEP (<25 GeV); QGSP (>12GeV) 
BIC (0-5 GeV); FTFB (>4 GeV) 
LEP (0-55 GeV); HEP (>25 GeV) 



Table 2. Indicates the physics model(s) which will be used for the primary pion interaction at each of the 
beam energies considered here. Where two models are indicated with a solidus, an appropriate random 
choice is made by GEANT 4. 



Physics List 


8 GeV 


10 GeV 


12 GeV 


15 GeV 


20 GeV 


30 GeV 


50 GeV 


80 GeV 


FTFP_BERT 


FTFP 


FTFP 


FTFP 


FTFP 


FTFP 


FTFP 


FTFP 


FTFP 


QGSP_BERT 


BERT 


LEP 


LEP 


LEP/QGSP 


LEP/QGSP 


QGSP 


QGSP 


QGSP 


QGSP_BERT_TRV 


BERT 


LEP 


LEP/QGSP 


QGSP 


QGSP 


QGSP 


QGSP 


QGSP 


QGSP_FTFP_BERT 


FTFP 


FTFP 


FTFP 


FTFP /QGSP 


FTFP/QGSP 


QGSP 


QGSP 


QGSP 


QGS_BIC 


LEP 


LEP 


LEP 


LEP/QGSB 


LEP/QGSB 


QGSB 


QGSB 


QGSB 


QGSP_BIC 


LEP 


LEP 


LEP 


LEP/QGSP 


LEP/QGSP 


QGSP 


QGSP 


QGSP 


FTF_BIC 


FTFB 


FTFB 


FTFB 


FTFB 


FTFB 


FTFB 


FTFB 


FTFB 


LHEP 


LEP 


LEP 


LEP 


LEP 


LEP 


LEP /HEP 


LEP /HEP 


HEP 



GEANT 4 also provides several interesting physics lists based wholly or in part on the CHIPS 
model [16]. In this picture, the result of a hadronic or nuclear interaction is the creation of a quas- 
mon (essentially an intermediate state of excited hadronic matter) which can dissipate energy by 
radiating particles statistically or by quark exchange with surrounding nucleons or clusters of nu- 
cleons. We have made studies of several CHI PS-based physics lists, but we choose not to present 
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Figure 3. Showing the principle of the cuts used to remove muons and non-interacting hadrons. The energies 
in the three calorimeters (in MIPs) are plotted for a 20 GeV data run. The three stacks in the ECAL are 
combined using weights proportional to the corresponding thicknesses of tungsten, i.e. 1:2:3. 



comparisons at this stage because, as discussed in Sect. 5.5, the model is still under active devel- 
opment. 

5. Selection of hadronic showers 

The principal task in the selection of events is to remove muons, electrons and protons from the 
sample. Muons are characterised by a small energy deposition in all three calorimeters. We show 
in figure || a scatter plot of the energies seen in the three calorimeters. A clear cluster close to 
the origin can be ascribed to muons. Bands are also seen corresponding to sharing of energy 
between ECAL and HCAL, and between HCAL and TCMT, with a few showers only starting in 
the TCMT. Accordingly, with a 10 GeV beam, we reject events for which the ECAL energy is less 
than 300 MIPs, the HCAL energy is less than 100 MIPs and the TCMT energy is less than 50 MIPs. 
These cuts vary linearly with energy, so that with a 30 GeV beam, for example, they are 329, 1 14 
and 64 MIPs respectively. Pions which do not interact in the ECAL are retained in the sample, to 
permit comparison of the interaction cross-section between data and simulation. The fraction of 
events removed as muons varies with beam energy, but it is typically ~ 5%. 

Double beam particles are very infrequent with the running conditions used in 2007, but in 
order to safely reduce any possible contribution from double beam particles, a cut is imposed, 
requiring that the total energy recorded in the ECAL+HCAL be less than 1 .5 times the beam energy. 
Fewer than 0.5% of events fail this cut. In order to reduce the influence of interactions upstream 




Figure 4. Energy observed in the ECAL for events with and without a signal from the beam Cerenkov 
counter. For this plot, the observed energy was converted from MIPS to GeV using a nominal conversion 
factor of 240 MIPs/GeV. (a) In the case of the 12 GeV negatively charged beam, the open region shows the 
total sample, and the shaded region the contribution with no Cerenkov signal, (b) For the 30 GeV positively 
charged beam the shaded histogram shows triggers with a Cerenkov signal, and the open ( red) histogram 
shows those with no Cerenkov signal. 



of the calorimeters, we remove events in which more than 50 MIPs are recorded in either of the 
first two layers of the ECAL. This cut removes ~ 2.5% of data, with only a weak dependence on 
beam energy. A cut against beam halo is also applied, requiring the shower centroid to lie within 
±50 mm in both x and y, corresponding approximately to the trigger acceptance; this cut also 
serves to remove a small number of showers close to the calorimeter edge, which are especially 
prone to lateral leakage. This cut removes ~ 5% of data at the lowest energy and ~ 1% at the 
highest energy, where the beam is narrowest. All of the cuts described so far are applied equally to 
data and Monte Carlo events. 

Electrons can provide a significant background to the pions in the negatively charged beams, 
and likewise protons are present in the higher energy positively charged beams. Their effect may 
be reduced by use of the Cerenkov counter. In figure ||(a) we compare, for typical negatively 
and positively charged beam runs, the distribution of the energy in the ECAL for triggers with and 
without a Cerenkov signal. In the case of the negatively charged beam we see a clear peak at the full 
beam energy in the events with a Cerenkov signal, which can be ascribed to electrons depositing 
their full energy in the ECAL. Therefore, for the negatively charged beams we demand that the 
Cerenkov have a signal. Conversely, for the positively charged beams used here, the Cerenkov 
pressure was set so that pions (and electrons) would yield a signal, while the heavier species would 
not. The ECAL energy distributions with and without a Cerenkov signal are shown in figure f|(b); 
we see that for the events with no Cerenkov signal, mainly protons, the distribution is shifted to 
lower values, and accordingly for these runs we demand that there be no Cerenkov signal. We also 
note that the contribution of positrons, which would show up as a peak close to 30 GeV in the 
sample of events with a Cerenkov signal, is small 3 . 

In our previous study of electron showers in the ECAL [j7|], it was necessary to remove showers 

3 From fitting the region around the full beam energy, we estimate a positron contribution of ~1% at 30 GeV, de- 
creasing to 0.4% at 50 GeV and 0.1% at 80 GeV. This has a negligible impact on the results shown in the present paper. 
No perceptible electron contribution is seen in the negatively charged beam samples. 
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Figure 5. Mean fraction of non-interacting pions (energy <100 MIPs) in the ECAL, plotted as a function 
of beam energy. For this purpose, the three ECAL stacks were combined with weights 1:1:1. The data 
are compared with the predictions of simulations using different GEANT4 physics lists. The models are 
separated into two plots in the interests of clarity, with the physics lists which incorporate the Bertini model 
on the left, and the others on the right. 



close to the edges of wafers, for two reasons. Firstly, the guard rings lead to a significant dead 
region (2 mm wide), which can lead to a significant fraction of the electron energy not being 
recorded. Secondly, we observed an effect of correlated crosstalk between the charge deposited in 
the guard rings and the peripheral cells of a wafer, whereby a square pattern of hits was sometimes 
seen in the core of a high energy shower [||]. The first of these problems is less important for 
hadronic showers because of their greater transverse width, and in any case the effect is sufficiently 
well simulated. The second problem is not seen in hadronic showers, because the energy densities 
in the core of the showers are much smaller. Therefore no further cuts on shower position are 
imposed. 



6. Comparison between data and simulation 
6.1 Total ECAL energy 

The total energy recorded in the ECAL is a useful starting point for comparison between data and 
simulation. The event selection cuts were designed to retain pion events wherever they interacted 
in the calorimeter system. A significant fraction of pions should not start to shower in the ECAL. 
These events are characterised by a MlP-like energy in all layers (apart from occasional 5 -ray 
emission) and accordingly the ECAL energy shows a large peak at ~50 MIPs. The fraction of 
such events can be used to test the interaction cross-sections in GEANT4. In figure || we show the 
fraction of selected pions depositing fewer than 100 MIPs in the ECAL, as a function of energy, 
compared with the simulations. Most of the models give a good description of the fraction of non- 
interacting pions at all energies, agreeing with data within 0.01-0.02. The LHEP physics list is the 
most discrepant. This gives confidence in the cross-sections simulated in GEANT4. It could also 
be interpreted as an indication that any residual beam contamination by kaons or (anti-)protons is 
small, since these species would be expected to have different interaction cross-sections. 

We now consider the energy deposited in the ECAL by those pions which have their first 
interaction in the ECAL. In the study of electromagnetic showers [f7|], it was appropriate to combine 
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Figure 6. Distributions of total energy recorded in the ECAL at 8, 15, 30 and 80 GeV (points with error 
bars), compared with Monte Carlo predictions using the QGSP_BERT physics list (solid histograms). The 
distributions are normalised to the same numbers of selected events (including the non-interacting peak). 



the energies in the three stacks of the ECAL taking account of their relative sampling fractions. 
Since the thicknesses of tungsten in the three stacks lie in the ratio 1:2:3, the weighting of the 
energy should be in roughly these ratios. It is less obvious that this is the correct procedure in the 
case of hadronic showers, but for the purposes of the present comparison, we choose to combine 
energies of each stack using the naive weighting factors of 1, 2 and 3. The level of agreement 
between data and simulation is not found to be sensitive to this choice. In figure ^ we compare the 
distribution of recorded energy in data with simulation using the QGSP_BERT physics list, at four 
typical energies. The non-interacting peak at low energies has been suppressed. The broad peaks 
in the distributions in figure ^ represent pions which started to shower in the ECAL, extending 
roughly to a point corresponding to all the energy being deposited in the ECAL (approximately 
250 MIPs per GeV). The main peak is quite well modelled at 8 and 15 GeV, while at 30 and 
especially 80 GeV the Monte Carlo predicts significantly more energy than observed in the data. 4 
In order to incorporate all models and energies into the analysis, a useful global measure is 
the average energy deposited in the ECAL, i.e. the mean values of the histograms such as those 
shown in figure |6[ excluding the non-interacting pions. This can be computed for data and for 
all the GEANT 4 physics lists under consideration, and is plotted in figure in the form of ratios 
of simulation to data. At 8 GeV, all of the models lie within 10% of the data, and most within 

4 In ref. ^ the energy response of the ECAL for electrons was seen to be linear up to 45 GeV. Since the energy 
densities in electromagnetic showers are significantly greater than those encountered in hadronic showers, we can be 
confident that saturation effects will not play any significant part in the hadron beam data. 
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Figure 7. Ratio of simulation to data for the mean energy recorded in the ECAL, plotted as a function of 
beam energy. The data are compared with the predictions of simulations using different GEANT 4 physics 
lists. 



5%. As energy increases, the LHEP model falls steadily further below the data. The other models 
either remain ~ 5 — 10% above the data or increase so that at high energies, all of the models 
lie ~ 5 — 10% above the data. Overall, FTF_BIC is the most consistent with this aspect of the 
experimental data. 

6.2 Transverse energy profile 

Good modelling of the transverse shower width is of importance for the development of particle 
flow algorithms, since it affects the degree of overlap between showers, and therefore the efficiency 
for separating them. For each hit in the ECAL, we determine the transverse distance between the 
centre of the pad and the shower barycentre in x and y. By histogramming this radial distance, we 
form the transverse shower profile. We weight the hits by their energy 5 , to emphasise the flow of 
energy in the shower. This also has the benefit of minimising possible residual effects of noisy 
cells. 

In figure |8] we show these transverse energy profiles for four typical energies, compared with 
the simulation based on the QGSP_BERT physics list. The entries in these plots are weighted by 
energy, and we want to focus here on the shape rather than the normalisation, which is essentially 
governed by the differences in energy response which were already considered in figure [7|. There- 
fore in figure || the distributions are normalised to unity, in order to facilitate the comparison of 
their shapes. We see that the data at 8 GeV are quite well modelled, while at higher energies the 
data tend to lie below the simulation at small radii, and above at higher radii; in other words the 
simulation underestimates the width of that part of the shower which is contained in the ECAL. 

In order conveniently to compare all models and energies, in figure || we show the mean 
energy-weighted shower radius (i.e. the mean values of distributions such as figure ||) as a function 
of energy. The observed shower becomes narrower with increasing energy, both in data and in all of 
the models. All of the physics lists underestimate the shower width at almost all energies, typically 
by around 10%. The models FTF_BIC and especially FTFP_BERT tend to lie closest to the data, 
especially at intermediate energies. 6 

5 In this case we do not apply different weighting factors to each stack. 

6 For comparison, the mean radius of electron-induced showers, selected using the methods of ref. M\, shows agree- 
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Figure 8. Radial distribution of hits (energy weighted) for data at four typical energies (points with errors) 
compared with Monte Carlo (solid histograms) using the QGSP_BERT physics list. The distributions are 
normalised to unity. 
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Figure 9. Mean energy-weighted shower radius in the ECAL as a function of beam energy. The data are 
compared with the predictions of simulations using different GEANT4 physics lists. 



Of course, the mean shower radius is only one measure of the transverse shower profile. In 
figure 10 we focus on the tails of the showers by plotting the radii needed to contain 90% or 
95% of the observed ECAL shower energy. As before, we find that most of the models tend to 
underestimate the data, and none of them gives a really satisfactory description. Again, the most 
successful physics lists are clearly FTFP_BERT and FTF_BIC. It should be stressed that these 
observations refer only to that part of the hadronic shower which is detected in the ECAL, i.e. 



ment between data and simulation at the ~ 1 — 2% or 0.1-0.2 mm level. 
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roughly the first interaction length — the transverse tails are not fully sampled, and nor is most of 
the longitudinal tail of the shower. 

6.3 Longitudinal distribution of interaction point 

The fine granularity of the CALICE ECAL makes it possible to identify the point at which the 
incident pion makes its primary interaction, and this will be used in some of the results shown 
below. In general one expects to see a track of MlP-like hits until the interaction point, followed 
by some multiparticle shower structure thereafter. The algorithm used for the present analysis has 
the merit of simplicity. Firstly, the energy in each layer is computed, after excluding isolated hits 
(those with no neighbours in adjacent cells in the same layer). The interaction layer is defined as 
the first layer containing at least 10 MIPs of energy, provided that at least two of the following 
three layers also show energy greater than 10 MIPs. In simulated events, the true interaction point 
is stored using information about the true activity in the Monte Carlo. In figure [n] we show the 
correlation between the true interaction layer and that found by our algorithm, for a typical energy 
and Monte Carlo physics list. The correct layer is identified within ±1 in ~70% and within ±2 
in almost 90% of events. The mean of the difference between the true and reconstructed layer lies 
within ±1 layer for all physics lists and energies studied here. 



In figure |12| we show the distribution of the reconstructed interaction layer for 30 GeV data 
compared with a typical Monte Carlo run. The alternation between odd and even layers reflects the 
additional material associated with PCBs and the carbon-epoxy support structure; passing from an 
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True Interaction Layer Reconstructed True Layer 

Figure 11. Left: comparison between reconstruction and truth for the layer identified as the interaction 
layer by the algorithm described in the text. This example corresponds to a 20 GeV 7t~ beam simulation, 
using the QGSP_BERT physics list. The grey scale indicates the number of events in each bin, where bins 
containing fewer than 10 events (amongst a total of more than 23000) have been suppressed for the sake of 
clarity. Right: distribution of (reconstructed — true) layer for this sample. 




Figure 12. Distribution of the reconstructed interaction layer in the ECAL for 30 GeV data (points), com- 
pared with Monte Carlo predictions using the QGSP _BERT physics list (solid histogram). 



odd-numbered layer to an even layer, a particle will encounter just a tungsten plate, while going 
from an even to an odd layer it will also pass through two PCBs and part of the mechanical support 
structure, adding an additional ~ 0.008Ai nt .. The attenuation of the beam through each stack is 
clearly seen, becoming stronger in later stacks as the tungsten thickness increases. The upward 
steps at layers 10 and 20 result from the increase of thickness of the tungsten plates at these points. 
The lower values in layers 1 and 2 are a by-product of the cut designed to remove upstream interac- 
tions. These features are all seen in the simulation, and are generally well modelled. The following 
analysis concentrates on interactions in the first stack, for which the agreement is particularly good. 
This distribution is similarly well described at all energies by all physics lists, which suggests that 
there is no significant problem with the high energy cross-sections on tungsten in this version of 
GEANT4. This is therefore not a directly useful distribution for discriminating between models. 
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Layer no. 



Figure 13. Energy per layer in the ECAL for 10 GeV electron data (points), compared with Monte Carlo 
predictions (solid histogram). 



6.4 Longitudinal energy profile 

It is of interest to study the longitudinal development of the showers. In the case of electrons, we 
find that the mean depth of electron-induced showers, selected using the methods of ref. [0], shows 
good agreement between data and simulation, to better than O.lXo (corresponding to less than 25% 
of a 1.4 mm tungsten layer), and the standard deviation is modelled to better than 0.03Xq. A typical 



comparison is shown in figure |13|. We therefore believe that the material content of the ECAL 
is well modelled. The case of hadronic showers is, however, complicated by the different points 
in the ECAL at which the shower is initiated. The observed longitudinal distribution of energy 
in the ECAL is the convolution of the intrinsic shape of showers with the distribution of shower 
starting points. In order to have some sensitivity to the composition of the shower in terms of 
different particle species, we would like to measure the shower profile with respect to the primary 
interaction point. 



In Sect. 5.3 we explained our algorithm for identifying the layer closest to the interaction 
point, and showed that it was reasonably reliable and well modelled. We now proceed to compute 
the shower profile in terms of layers after the interaction layer. This is made more complicated by 
the different sampling fractions in the three ECAL stacks. We circumvent this by the following 
expedient. For hits in the first stack (tungsten thickness 1.4 mm), we simply use the measured 
energy in each layer. In the second stack (tungsten thickness 2.8 mm), in addition to the energy 
measured in each layer, we introduce an additional fictitious pseudolayer midway between each 
physical layer, whose energy is estimated by linear interpolation between the layers on either side. 
Similarly, in the third stack (tungsten thickness 4.2 mm) we introduce two pseudolayers between 
each physical layer. In this way, the energy deposition is estimated in 60 layers each separated 
by an effective separation of 1.4 mm of tungsten. The longitudinal shower profile (i.e. the mean 
energy, in MIPs, per layer) is then trivially computed starting at the interaction layer measuring 
the depth in units of these effective 1.4 mm layers. This calculation will lead to some non-trivial 
correlations between neighbouring bins in the longitudinal profile. It also neglects differences in 
the material between successive samplings caused by other passive materials such as PCBs, though 
this effect will be averaged out by combining different interaction layers. Furthermore, both these 
features will also be present in the simulation. For this study, we restrict the interaction layer to lie 
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in the first stack, so that at least 50 effective layers are available for observation the shower profile, 
corresponding to ~ 20 Xo or 0.8 Aj nt .. 

In the case of the GEANT4/Mokka simulation, it is possible to store a breakdown of the frac- 
tions of the energy of each hit contributed by various particle species. In this way, one can break 
down the longitudinal shower profile into its constituent parts. We choose to separate the energies 
deposited by electrons, positrons, "mesons" (charged pions, kaons and muons), protons (including 
antiprotons) and "others" (e.g. deuterons, a -particles, hyperons). 

In figure |l4| we show, for 12 GeV n~ , the longitudinal energy profiles data compared with 
simulations using different physics lists. In each case, we also show the simulation broken down 
into each of the above categories. The "mesons" show a smooth slowly falling shape, while the 
electrons and positrons show the characteristic rise and fall of an electromagnetic shower, peaking 
after ~10-15 layers (~ 4 — 6Xo). The electrons, however, also exhibit a long tail; this is because 
part of the ionisation energy loss of other particles is simulated in the form of discrete knock- 
on electrons, and therefore is included in the electron contribution. The protons, and the small 
contribution from "others", seem to exhibit two components - a sharp peak in the first few layers 
which, because of its short range, can be attributed to low energy nuclear break-up fragments, and 
a long range tail similar to the mesons. We note that there are significant differences between 
the models, most conspicuously in the short range proton component. The comparison with data 
demonstrates that our calorimeter has sufficient longitudinal granularity to be able to offer some 
discrimination between these components. 



The main features of the data distribution in figure 14 may be summarised as follows: 



• A small peak is seen in the first few layers, which we attribute to nuclear fragments, mainly 
protons. As noted above, the models exhibit considerable variation in their predictions in this 
region, associated with differences in their proton yields. None agrees well with data, with 
the FTF-based physics lists overestimating the data, and the others undershooting. 

• A broad peak occurs around layer 10, which we ascribe to the electromagnetic component. 
This is tolerably well modelled by all physics lists, with some variations in normalisation. 

• A long tail follows, which appears to be generally well modelled. 

In figure [l| we show similar distributions for two physics lists, FTFP_BERT and QGSP_BERT, 
at four typical energies. Clearly both models have imperfections, but on balance FTFP_BERT is 
probably the more successful. 

In order to quantify these observations, and extend them to all energies and models, we show 
in figure |l(] ratios of simulation to data as a function of energy for three regions of the longitudinal 
profile: layers 1-3 (where the contribution from nuclear breakup dominates), layers 5-20 (domi- 
nated by the showering of electrons and photons) and layers 30-50 (the tail, dominated by pene- 
trating hadrons). This makes it clear that the greatest differences between the models occur in the 
first few layers, where discrepancies up to ~ 20 % are observed. The FTFP_BERT and FTF_BIC 
models lie consistently above the data, while LHEP, QGSP_BIC and QGS_BIC he consistently 
below. The other physics lists make a transition between the two regimes, as their model content 
changes with energy. In the electromagnetic-dominated region, the most obvious outlier is LHEP; 
the other models all agree with data within ~ 15%, with the FTFP_BERT and FTF_BIC models 
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Figure 14. Longitudinal energy profiles for 12 GeV %~ data {shown as points), compared with simulations 
using different physics lists. The mean energy in MIPs is plotted against the depth after the initial interaction, 
in units of effective 1.4 mm tungsten layers. The total depth shown corresponds to ~ 20 Xq or 0.8 Aj nt .. The 
breakdown of the Monte Carlo into the energy deposited by different particle categories is also indicated. 




Figure 15. Longitudinal energy profiles for data (shown as points) compared with simulations using two 
physics lists, QGSP_BERT and FTFP_BERT, at four typical energies. The breakdown of the Monte Carlo 
into the energy deposited by different particle categories is also indicated. 



giving the best description. In the tails, most models lie within ~10% of data; LHEP is consistently 
low, as is FTF_BIC at lower energies. 

On balance, it appears that the FTFP_BERT physics list, while not perfect, gives the best 
overall description of the longitudinal development of these showers. We emphasise, however, that 
this remark refers only to the early part of the shower which is developed in the ECAL; we are not 
sensitive to the later parts of the shower. 
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6.5 Physics lists under development 

There are several physics lists in GEANT 4.9.3 which use the CHIPS model JI^ ] alone or in con- 
junction with other models. This model has undergone substantial development recently, and is 
now capable of modelling all parts of the hadronic interaction process. We made a number of 
studies of these physics lists using a j3-test version 4.9.3.b01 of GEANT. The results were very en- 
couraging, and one of the physics lists tested, QGSC_CHIPS, gave a better description of our data 
than any of the physics lists discussed here. However, the CHIPS model is still being tuned, and in 
the released version of GEANT 4.9.3, the CHIPS-based models are less successful when confronted 
with our data. It would be premature to show results and draw conclusions while development is 
ongoing, but the CHIPS-based models seem to be an interesting and promising new avenue, and it 
is clear that our data have the power to discriminate between tunings of such models. 

7. Summary 

We have studied showers induced by charged pions in the CALICE silicon-tungsten ECAL. The 
calorimeter has high transverse and longitudinal granularity, which allows us to study the early part 
of the shower development in unprecedented detail. Using pions in the energy range 8-80 GeV, 
we have compared the data with GEANT 4 simulations using eight different physics lists. Several 
observables were examined - the total energy deposited in the ECAL, the distribution of interac- 
tion points, the transverse shower energy profile, and the longitudinal profile with respect to the 
interaction point. The most sensitive observables seem to be the shower profiles. The data tend to 
show a greater shower width than the simulation for most energies and physics lists. The closest 
description of the transverse profiles is provided by the FTFP_BERT and FTF_BIC physics lists. 

The longitudinal profile seems the most interesting observable, because the three main compo- 
nents of the shower induced by the primary interaction (energetic hadrons, photons and low energy 
nuclear fragments) can be, to some extent, distinguished through their different rates of shower de- 
velopment. An ECAL using tungsten is particularly useful in this regard, because of its small ratio 
of Xq I Aint. , which amplifies the differences in shower development between various components of 
the shower. It would be naive to expect any of the physics lists to give a perfect description of the 
data, but it seems clear from our study that the LHEP physics list has serious deficiencies, and that, 
in the framework of the current version 4.9.3 of GEANT 4, FTFP_BERT list is the most successful. 
It is to be hoped that these observations will provide useful guidance to the GEANT 4 developers. 
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Figure 16. Ratio of simulation to data for three different regions of the longitudinal energy profile: ( top pair 
of plots) layers 1-3, dominated by nuclear breakup; ( centre pair) layers 5-20, dominated by electromagnetic 
showers; and (bottom pair) layers 30-50, dominated by penetrating hadrons. 
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